
% adjusted wage=w*E^(1/epsilon) 

clc
clear all

w=ones(26,2); % no wage differentials
load moments_from_Seoul_survey.mat
load distance_GIS.mat

% w: 25 district + home sector
H_R(:,1)=H_R_young; % young
H_R(:,2)=H_R_old;   % old

H_M_target(:,1)=[H_M_target_young_weekday*(1-0.17); 0.17*sum(H_R_young)]; % young weekdays
H_M_target(:,2)=[H_M_target_young_weekend*(1-0.46); 0.46*sum(H_R_young)]; % young weekends
H_M_target(:,3)=[H_M_target_old_weekday*(1-0.61); 0.61*sum(H_R_old)];   % old weekdays
H_M_target(:,4)=[H_M_target_old_weekend*(1-0.67); 0.67*sum(H_R_old)];   % old weekends

nu(1,1)=0.1413; % weekdays
nu(1,2)=0.1666; % weekends
ncity=25;

epsilon(1,1)=4.1642; % weekdays
epsilon(1,2)=4.9144; % weekends

global w H_R H_M_target nu d_GIS ncity

% Tolerance;
eps = 1e-3;

%% transformed wage (omega)
guess = ones(26,1);
LBD=10^(-10)*ones(26,1);
options = psoptimset('Display','iter','TolFun',eps,'TolX',eps);
[x1] = patternsearch(@Calibrate_omega_function_young_weekday,guess,[],[],[],[],LBD,[],[],options);
omega(:,1)=x1;
guess = ones(26,1);
LBD=10^(-10)*ones(26,1);
options = psoptimset('Display','iter','TolFun',eps,'TolX',eps);
[x2] = patternsearch(@Calibrate_omega_function_young_weekend,guess,[],[],[],[],LBD,[],[],options);
omega(:,2)=x2;

guess = ones(26,1);
LBD=10^(-10)*ones(26,1);
options = psoptimset('Display','iter','TolFun',eps,'TolX',eps);
[x3] = patternsearch(@Calibrate_omega_function_old_weekday,guess,[],[],[],[],LBD,[],[],options);
omega(:,3)=x3;
guess = ones(26,1);
LBD=10^(-10)*ones(26,1);
options = psoptimset('Display','iter','TolFun',eps,'TolX',eps);
[x4] = patternsearch(@Calibrate_omega_function_old_weekend,guess,[],[],[],[],LBD,[],[],options);
omega(:,4)=x4;

%% E
E_weekday(:,1)=omega(:,1)./(w(:,1).^epsilon(1,1));
E_weekday(:,1)=E_weekday(:,1)/geomean(E_weekday(:,1)); % geometric mean 1
E_weekend(:,1)=omega(:,2)./(ones(26,1).^epsilon(1,2));
E_weekend(:,1)=E_weekend(:,1)/geomean(E_weekend(:,1)); % geometric mean 1
E_weekday(:,2)=omega(:,3)./(w(:,2).^epsilon(1,1));
E_weekday(:,2)=E_weekday(:,2)/geomean(E_weekday(:,2)); % geometric mean 1
E_weekend(:,2)=omega(:,4)./(ones(26,1).^epsilon(1,2));
E_weekend(:,2)=E_weekend(:,2)/geomean(E_weekend(:,2)); % geometric mean 1

save ('Output/calibrated_E.mat','E_weekday','E_weekend')
